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ABSTRACT 

We study the impact of baryons on the distribution of dark matter in a Milky Way-size halo by 
comparing a high-resolution, moving-mesh cosmological simulation with its dark matter-only 
counterpart. We identify three main processes related to baryons - adiabatic contraction, tidal 
disruption and reionization - which jointly shape the dark matter distribution in both the main 
halo and its subhalos. The relative effect of each baryonic process depends strongly on the 
subhalo mass. For massive subhalos with maximum circular velocity Umax > 35 km s“^, adi¬ 
abatic contraction increases the dark matter concentration, making these halos less susceptible 
to tidal disruption. For low-mass subhalos with Umax < 20 kms“^, reionization effectively 
reduces their mass on average by « 30% and Umax by « 20%. For intermediate subhalos 
with 20 kms“^ < Umax < 35 kms“^, which share a similar mass range as the classical 
dwarf spheroidals, strong tidal truncation induced by the main galaxy reduces their Umax- As 
a combined result of reionization and increased tidal disruption, the total number of low-mass 
subhalos in the hydrodynamic simulation is nearly halved compared to that of the Wbody 
simulation. We do not find dark matter cores in dwarf galaxies, unlike previous studies that 
employed bursty feedback-driven outflows. The substantial impact of baryons on the abun¬ 
dance and internal structure of subhalos suggests that galaxy formation and evolution models 
based on Wbody simulations should include these physical processes as major components. 
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1 INTRODUCTION 

A major achievement in observational cosmology is the discovery 
that our Universe is composed of ~ 4% baryons, 20% dark matter 
(DM), and 76% dark energy (DE) (Frieman et al. 2008). The first 
observational evidence for DM dates back to 1933 when Zwicky 
noted a missing mass problem in the Coma cluster of galaxies: the 
visible galaxies account for only a small fraction of the total mass 
inferred from the dynamics (Zwicky 1937). More evidence came 
later from galactic rotation curves in spiral galaxies (Rubin et al. 
1980), gravitational lensing, and the Bullet cluster which shows 
an offset of the center of the total mass from that of the baryons 
(Clowe et al. 2006). The first compelling observational evidence for 
DE was found later in 1998 when two teams studying Type la su¬ 
pernovae independently found that the expansion of the universe is 
accelerating (Riess et al. 1998; Perlmutter et al. 1999). This finding 
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has been confirmed by subsequent supernovae observations, and 
independent evidence from galaxy clusters (e.g., Vikhlinin et al. 
2006; Allen et al. 2011), large-scale structure (e.g., Tegmark et al. 
2006; Addison et al. 2013), and the cosmic microwave background 
(e.g., Spergel et al. 2007; Komatsu et al. 2011; Hinshaw et al. 2013; 
Planck Collaboration et al. 2015). 

These observations motivate the current “standard model” 
of cosmology (ACDM), where dark energy and cold dark mat¬ 
ter shape the formation and evolution of cosmic structures (e.g., 
Frenk & White 2012; Kravtsov & Borgani 2012; Conselice 2014; 
Somerville & Dave 2014). To date, numerous ACDM cosmologi¬ 
cal simulations have produced clumpy and filamentary large-scale 
structures as seen in galaxy surveys (e.g., Navarro et al. 1997; 
Springel et al. 2005b, 2006; Gao et al. 2012; Vogelsberger et al. 
2014a), and have confirmed that structures form through hierar¬ 
chical assembly in CDM dominated universes. However, on small 
scales (i.e. less than ~ 10 kpc), there appear to be a number of 
tensions between predictions from the ACDM model and observa- 
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tions, notably: (1) the “missing satellites problem”, in which the 
abundance of subhalos produced by A^-body simulations is orders 
of magnitude larger than the two dozens of satellites observed in 
the MW (e.g., Klypin et al. 1999; Moore et al. 1999; Kravtsov et al. 
2004; Kravtsov 2010); (2) the “too big to fail problem”, in which 
Al-body simulations produce overly dense massive subhalos com¬ 
pared to the brightest dwarf galaxies in the MW and Local Group 
(Boylan-Kolchin et al. 2011, 2012; Tollerud et al. 2014; Garrison- 
Kimmel et al. 2014a); and (3) the “core-vs-cusp problem”, in which 
the central dark matter density profiles of DM-dominated dwarf 
spheroids are observed to apparently feature smooth cores instead 
of the cusps that are generically predicted by CDM models (e.g., 
Gilmore et al. 2007; Evans et al. 2009; de Blok 2010; Amorisco 
& Evans 2012; Strigari et al. 2010; Martinez 2013). These prob¬ 
lems have motivated alternative models such as self-interacting DM 
(Dave et al. 2001; Vogelsberger et al. 2012; Elbert et al. 2014; Vo- 
gelsberger et al. 2014b), or warm DM (e.g., Polisensky & Ricotti 
2014; Schneider et al. 2014; Kennedy et al. 2014). 

However, some of the discrepancies were reported in DM-only 
simulations in which the dynamical coupling between baryons and 
dark matter was ignored. While this assumption may be justified 
on large scales, this is no longer the case on kpc scales, where the 
density starts to be dominated by baryons, and the dynamics be¬ 
comes governed by baryonic processes such as gas dynamics, star 
formation, black hole accretion, and feedback from stars and active 
galactic nuclei (AGN). 

An impact of baryons on the dark matter arises from different 
spatial distributions of the two components. A well-known effect 
is adiabatic contraction (Young 1980; Barnes & White 1984; Blu- 
menthal et al. 1986; Ryden & Gunn 1987; Gnedin et al. 2004, 2011; 
Zemp et al. 2012; Pillepich et al. 2014), which causes an increase of 
the mass concentration of dark matter in the center of a galaxy due 
to gas inflow as a result of cooling. The increased DM mass concen¬ 
tration and the presence of a stellar disk can produce stronger tidal 
forces, which have been suggested as a way to significantly affect 
the abundance and distribution of subhalos and satellite galaxies in 
the Milky Way (MW) (e.g. D’Onghia et al. 2010a; Pefiarrubia et al. 
2010; Zolotov et al. 2012; Arraki et al. 2014). 

Hydrodynamic simulations have become powerful tools to in¬ 
vestigate the response of dark matter to baryons and vice-versa, 
thanks to recent progress in numerical methods (e.g. Springel 2010; 
Read & Hayfield 2012; Hopkins 2013; Hu et al. 2014) and physi¬ 
cal modeling (Vogelsberger et al. 2013; Stinson et al. 2013; Aumer 
et al. 2013; Hopkins et al. 2014) that improved upon long-standing 
issues in the field (Agertz et al. 2007; Sijacki et al. 2012; Torrey 
et al. 2012). Recent cosmological simulations such as the Illustris 
(Vogelsberger et al. 2014a) and EAGLE (Schaye et al. 2015) projects 
were able to reproduce different galaxy populations that resemble 
the observed ones both locally and in the high-redshift Universe 
(e.g. Vogelsberger et al. 2014c; Genel et al. 2014). In particular, 
Vogelsberger et al. (2014a) showed that subhalos in DM-only sim¬ 
ulation are more prone to tidal disruption than those in hydrody¬ 
namic simulations, leading to a depletion of satellites near galaxy 
cluster centers and a drop in the matter power spectrum on small 
scales. 

Several solutions have been proposed to solve the “too big to 
fail” problem, including tidal effects (Zolotov et al. 2012; Arraki 
et al. 2014), which however might be insufficient in some cases 
(for instance in M31, see Tollerud et al. 2014), a mass-dependent 
abundance of subhalos which may alleviate the problem if a lower 
total mass of the MW is assumed (Wang et al. 2012; Vera-Ciro 
et al. 2013; Sawala et al. 2014a), and strong outflows driven by 


supernovae explosions which can have a direct impact on the cen¬ 
tral dark matter content in dwarf galaxies, possibly leading to a 
cored profile (e.g. Navarro et al. 1996; Governato et al. 2012; 
Teyssier et al. 2013; Madau et al. 2014; Brooks & Zolotov 2014; 
Ogiya & Burkert 2015). However, Boylan-Kolchin et al. (2012) 
and Garrison-Kimmel et al. (2013) argue that the required energy 
from supernovae may not be sufficient given the low stellar mass 
in some of the dwarf galaxies (but see also the energy argument 
by Madau et al. 2014). Moreover, most hydrodynamic simulations 
(e.g. Mashchenko et al. 2008; Governato et al. 2012; Teyssier et al. 
2013; Madau et al. 2014) are focused on dwarf galaxies in field en¬ 
vironments, which may not be representative for the dwarfs in the 
MWorM31. 

In order to investigate baryonic effects on dark matter in dwarf 
galaxies of the MW, we need high-resolution, cosmological hydro- 
dynamic simulations which produce a spiral galaxy with properties 
similar to those of the MW. Producing MW-like disk galaxies in 
cosmological simulations has been a decade-long challenge, but re¬ 
cently several groups have succeeded in this endeavor (Agertz et al. 
2011; Guedes et al. 2011; Aumer et al. 2013; Okamoto 2013; Hop¬ 
kins et al. 2014). Equipped with the same implementation of baryon 
physics as in the Illustris simulations, Marinacci et al. (2014a) suc¬ 
cessfully produced MW-size disk galaxies in a suite of zoom-in 
simulations. These simulations used the same initial conditions 
as the Aquarius Project (Springel et al. 2008), and the highest- 
resolution hydrodynamical rerun (Aq-C-4) has sufficient resolution 
to identify and study the formation history and properties of the 
predicted dwarf galaxies. 

In this work, we use both DM-only and hydrodynamical sim¬ 
ulations of Aq-C-4 by Marinacci et al. (2014a) to study the impact 
of baryon processes on the halo/subhalo properties and the sub¬ 
halo abundance. We will not limit our study on the bright satellites 
alone, i.e. those subhalos containing stars, but we will also analyze 
the “dark” ones. As it turns out, even the “dark” subhalos are sys¬ 
tematically affected by baryonic processes in terms of their spatial 
distribution and mass functions. 

This paper is organized as follows. In Section 2, we describe 
the numerical technique used in our simulations and the structure 
identification. The impact of baryons on the smooth dark matter 
distribution in the main halo and the global statistics of subhalos 
are presented in Section 3. In Section 4, we investigate the impact 
of baryons on the total mass, the DM density profiles and Umax val¬ 
ues of objects extracted from a matched subhalo catalogue of the 
DM and Hydro simulations. We aim to determine the main physi¬ 
cal processes that shape the DM content in subhalos by tracking the 
assembly history and evolution of bright satellites and “dark” sub¬ 
halos. We discuss the implications of our study and its limitations 
in Section 5, and summarize our main findings in Section 6. 

2 METHODS 
2.1 The Simulations 

In this study, we use two cosmological simulations of a MW-size 
halo, one being the full hydrodynamical Aq-C-4 run by Marinacci 
et al. (2014a) (referred to as “Hydro” hereafter), and the other 
being a control DM-only simulation of the same halo (referred 
to as “DM0” hereafter). This Aq-C halo was selected as a close 
match to the MW for the Aquila Comparison Project (Scanna- 
pieco et al. 2012), as well as several other studies (Wadepuhl & 
Springel 2011; Okamoto 2013; Sawala et al. 2012). The hydrody- 


(© 2015 RAS, MNRAS 000, 1-22 



Baryonic impact on dark matter 3 


namical Aq-C-4 simulation by Marinacci et al. (2014a) was per¬ 
formed with the moving mesh code Arepo (Springel 2010). The 
simulation adopted a physical model for galaxy formation and evo¬ 
lution developed by Vogelsberger et al. (2013), which includes su¬ 
pernovae feedback, metal enrichment and stellar mass return, AGN 
feedback (Di Matteo et al. 2005; Springel et al. 2005a; Sijacki et al. 
2007), and a spatially uniform, redshift-dependent ionizing back¬ 
ground by Faucher-Giguere et al. (2009), which leads to complete 
re-ionization of neutral hydrogen by z = 6. Thermal feedback from 
supernovae was implemented following a hybrid ISM model devel¬ 
oped by Springel & Hernquist (2003), and galactic outflows were 
launched with a velocity scaled with the local dark matter veloc¬ 
ity dispersion of the host halo, following a kinetic model similar to 
Okamoto et al. (2010) and Puchwein & Springel (2013). 

The Aq-C-4 Hydro simulation has a mass resolution of 5.0 x 
10"* Mq for gas and stars, and 3.0 x 10® Mq for the DM compo¬ 
nent (see Table 1 in Marinacci et al. 2014a), sufficient to simul¬ 
taneously follow the main galaxy and its classical dwarf galaxies 
with a maximum circular velocity Wmax between 12 and 24 km 
(Boylan-Kolchin et al. 2011). The gravitational softening length in 
the high-resolution region was kept fixed in comoving coordinates, 
corresponding to a physical length of 340 pc at 2 = 0. We re-run 
a DM-only simulation of Aq-C-4 with the same numerical param¬ 
eters controlling the force and time integration accuracy as used in 
Marinacci et al. (2014a). Thus, the effects of baryonic processes on 
the dark matter distribution can be well studied by comparing the 
DM simulation and its Hydro counterpart. At 2 = 0, the proper¬ 
ties of the central galaxy of Aq-C-4 in the Hydro run are in very 
good agreement with those of a typical disk-dominated galaxy in 
terms of the mass budget in various components, the morphology, 
and the star formation history (Marinacci et al. 2014a). The proper¬ 
ties of the diffuse gas and the metal distribution are also consistent 
with observations (Marinacci et al. 2014b). Moreover, we note that 
the robustness of the results was verified by a resolution study in 
Marinacci et al. (2014a). 


2.2 Structure Identification 

To identify subhalos both in the DM0 and Hydro simulations, the 
snapshots were post-processed with the Amiga Halo Finder (AHF, 
Knollmann & Knebe 2009).^ The AHF algorithm identifies struc¬ 
tures based on density estimates calculated with an adaptive refine¬ 
ment technique, and naturally builds a halo-subhalo-subsubhalo hi¬ 
erarchy. The extent of a halo is determined by its density, p(rvir) = 
Avir(2)pbg, where p(rvir) is the mean density within the virial ra¬ 
dius Tvir, pbg is the background density, and Avir(2) = 178 is the 
adopted virial overdensity. 

AHF performs an iterative process to remove unbound parti¬ 
cles until the final result converges to a set of bound particles within 
Tvir. These sets of particles form the halo and subhalos. The code 
then calculates various properties of the halos and subhalos, such 
as the mass in different components, the maximum value of the ro¬ 
tation curve Wmax, and the spin parameter. The results of AHF and 
other substructure finders such as SUBFIND (Springel et al. 2001) 


^ The code is available at http://popia.ft.uam.es/AHF/ 
Download.html. In this study, we use the version ahf-vl.0-084. It also 
contains an analysis tool called MergerTree, which we have used to con¬ 
struct the merger tree and to cross-match the subhalos between the DM0 
and Hydro simulations. 


are generally in very good agreement (e.g. Onions et al. 2012; Pu¬ 
jol et al. 2014), and any residual differences are not expected to 
influence our results. 

The IDs of collisionless particles are preserved in our simula¬ 
tions since there is no mass exchange between them. This allows us 
to construct subhalo merger trees using the built-in module Merg¬ 
erTree of AHF, which relies on tracking the membership of dark 
matter and star particles (identified by their IDs) within the different 
halos and subhalos. The merger trees are constructed for both the 
DM and Hydro simulations. For each halo/subhalo, we only con¬ 
sider the most massive progenitor in the previous snapshot as its 
parent. In addition, we only consider the halos/subhalos comprised 
of at least 10 particles. We have carefully checked the validity of 
the constructed merger trees by visually comparing the evolution¬ 
ary paths of each individual object in terms of its position, velocity 
and mass. There are rare cases when a subhalo is not detected by 
AHF in one snapshot output when the subhalo closely passes the 
center of its host halo. These objects usually reappear in the next 
AHF catalogue if they have not been disrupted at pericenter. To 
avoid complications, we discard such subhalos in this study. 

With the MergerTree analysis package, we can also cross¬ 
match the 2 = 0 snapshots of the DM and Hydro simulations using 
the dark matter particles. We verify that this cross-match between 
the Hydro and DM0 simulations is able to identify the “same” ob¬ 
jects by comparing their evolutionary paths. In Section 4.1, we 
show the orbital and mass growth histories of several of these 
matched objects. We note however that substantial orbital phase 
offsets are expected to appear in most pairs due to the inclusion of 
baryonic processes; thus the positions of the subhalos in the two 
simulations are not expected to exactly match each other. Also, 
some halos/subhalos identified in the DMO simulation do not have 
counterparts in the Hydro simulation, given that substructures are 
destroyed at a higher rate in the latter run. 


3 BARYONIC IMPACT ON THE PROPERTIES OF THE 
MAIN GALAXY AND ITS SATELLITES 

One of the goals of this study is to identify the main physical pro¬ 
cesses shaping the distribution of dark matter in the main galaxy 
and its substructures. In this section, we focus on important galaxy 
properties such as the spatial distribution, abundance, and mass 
function of satellites, as well as the DM density profiles in the DMO 
and Hydro simulations. 

3.1 Spatial Distribution and Abundance of Subhalos 

In Figure I, we show projected dark matter density maps at 2 = 0 
for a slice of thickness 250 /i“^kpc centered on the MW-size halo 
for the DMO (top panel) and the Hydro (bottom panel) simula¬ 
tions, respectively. Numerous substructures are clearly visible in 
both simulations. Despite the overall similarity in the morphology 
and size of the main halo between the two simulations, there are 
notable differences in the abundance and spatial distribution of the 
subhalos, especially in the central region, as demonstrated in Fig¬ 
ure 2, which shows the positions of all the DM subhalos and bright 
satellites (subhalos that contain stars) within the virial radius of the 
main galaxy at the present epoch. It is clear that there are fewer sub¬ 
halos in the Hydro simulation than in the DMO one, in particular, 
the large number of low-mass subhalos found in the DMO simula¬ 
tion is clearly reduced in the Hydro case. Moreover, it is seen that 
only a fraction of the subhalos presented in the DMO simulation 
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Figure 1. Projected dark matter density maps within a 250 h~^kpc slice 
centered on the main halo at redshift z = 0 in the DM0 (top panel) and 
Hydro (bottom panel) simulations, respectively. The size of the displayed 
region is 0.7 h~^Mpc on a side. 




Figure 2. The spatial distribution of subhalos at redshift 2 = 0 in the DM0 
(top panel) and Hydro (bottom panel) simulations, respectively. The black 
filled circles represent dark matter subhalos, while the red filled circles rep¬ 
resent bright satellites which have formed stars. The size of the symbols is 
scaled with the subhalo mass. The solid yellow circle indicates the virial 
radius calculated by AHF (with overdensity Avir = 178). The Hydro sim¬ 
ulation produces fewer subhalos than the DM0 counterpart, with a pro¬ 
nounced depletion of low-mass subhalos near the central region. The bright 
satellites are only a small fraction of the entire subhalo population. 


can be found in the Hydro simulation. Note that not necessarily the 
most massive ones are able to host bright satellites that form stars. 

To quantify the spatial distribution of subhalos, we compare 
the radial number density of subhalos in different mass ranges from 
both the DM0 and Hydro simulations in Figure 3. It has been 
shown that in A-body simulations the spatial distribution of sub¬ 
halos follows a universal function which is less concentrated than 
the density profiles of dark matter halos. It can be parameterized by 
the Einasto profile (Gao et al. 2012): 


n{r)/{n)uMO = (n)_2,oexp 



( 1 ) 


We fit our data with the Einasto profile for subhalos in two 
mass ranges, delineated by maximum circular velocities Umax > 
5kms“'^ and Umax > lOkms”'^ as shown in Figure 3. The ra¬ 
dial abundance of subhalos in the Hydro simulation is consistently 
lower than that in the DM0 one, and the effect increases towards 
the center for the most massive subhalos. There are not enough data 
points within 0.2rvir for Umax > 20 km to perform reliable fit¬ 
ting with an Einasto profile. The cumulative radial distribution of 
subhalos in Figure 4 confirms this trend. The total number of subha¬ 
los at any given radius, with the exception of the innermost regions 
(u/uvir 0.1) in the Umax > 5kms“^ cut, is consistently lower in 
the Hydro simulation than in the DM0 run, and this is particularly 
evident for the most massive subhalos in the central regions. 
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Figure 3. The number density of subhalos in different mass ranges as a 
function of distance to the center of the main halo, both for the DMO 
(black symbols) and Hydro (red symbols) simulations. The distance r is 
normalized to the virial radius of the main halo, while the subhalo abun¬ 
dance is normalized to (n,)i3MO> the total number of DM subhalos identi¬ 
fied in the DMO simulation divided by the entire volume enclosed by rvi, ■ 
The top, middle and bottom panels show subhalos in three different mass 
ranges, as indicated by the maximum circular velocity Umax > 5 km 
Umax > 10 km s“t and Umax > 20 km s“t respectively. The error bars 
are computed using the Poisson error \/Nr, where Nr is the number of sub¬ 
halos within each radial bin. The dashed lines are fits to the Einasto profile, 
as given by Eqn. (1). 


Figure 4. The cumulative number of subhalos in different mass ranges as 
a function of distance to the center of the main halo, both for the DMO 
(black lines) and Hydro simulations (red lines). The three mass ranges are 
the same as in Figure 3: Umax > 5kms“^, Umax > 10kms“^, and 
Umax > 20kms“^, respectively. 


Both Figure 3 and Figure 4 suggest that subhalos are subject 
to being disrupted more easily in the Flydro simulation. A simi¬ 
lar radial distribution can also be found in D’Onghia et al. (2010a) 
and Yurin & Springel (2015), in which it was suggested that the 
reduction of subhalos was due to enhanced tidal effects and accel¬ 
erated disruption rates from a combination of DM contraction and 
the presence of the stellar disk. In addition, enhanced dynamical 
friction from the adiabatically contracted dark matter distribution 
of the main halo would cause the subhalos to sink more rapidly. 
The combination of these factors results in fewer massive subhalos 
in the central region in the Hydro simulation than in the DMO one. 
We will address the impact of the central galaxy on the abundance 
of its satellites in Section 4. 

A comparison of the cumulative distribution of DM subhalos 
between the DMO and Hydro simulations is shown in Figure 5. 
Both simulations show a power-law distribution of the subhalo 
abundance, N{> Umax) oc Umax in terms of maximum circular 
velocity, and N{> Msub) oc M~^ in terms of mass, similar to the 
relations reported for DM subhalos based on the Aquarius simula¬ 
tions (Springel et al. 2008) and the Phoenix simulations (Gao et al. 
2012). The slopes of both distribution functions are similar in our 
DMO and Hydro simulations. However, the total number of subha¬ 
los in the Hydro simulation is consistently lower by ~ 50% than in 
the DMO case, except for the range where Umax > 35 kms“^ (or 
equivalently Msub > 4 x 10® Mq in terms of mass). 

The bright satellites, which are here defined as subhalos con¬ 
taining stars in the Hydro simulation, show a different distribution 
from the DM subhalos at a critical point of Umax ~ 20 kms“^ 
(corresponding to a mass of Msub ~ 10® M©). At the low-mass 
end, the probability of a subhalo hosting stars steadily decreases as 
Umax decreases. The “missing satellite problem” appears clearly 
striking if we simply compare the number of DM subhalos at 
Umax < 10 kms“^ with observations of Penarrubia et al. (2008a), 
because the former is more than two orders of magnitude higher. 
However, the number of satellites (i.e. subhalos with stars) is much 
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Figure 5. . The cumulative distribution of the number of subhalos from the 
DM0 (black solid line) and Hydro (red solid line) simulations, as a function 
of the maximum circular velocity i^max (top panel) and the subhalo mass 
Mgub (bottom panel). The gray dashed lines are fits from the literature, 
N{> Umax) OC UmL (top panel), or N{> M^ub) « (bottom 

panel). Bright satellites (subhalos that have stars) are represented by the 
pink solid curve, while observations by Penarrubia et al. (2008a) are shown 
with blue dots, for comparison. 

closer to the observations, and the discrepancy between the two 
becomes even smaller when detection and completeness limits of 
current surveys are accounted for. 

At the massive end, Umax > 20 kms“^, the number of bright 
satellites agrees well with observations and it matches that of DM 
subhalos. The value of Umax ~ 20 kms“^ marks a transition in 
dwarf galaxy formation shaped by reionization, similar to previ¬ 
ous studies (Okamoto et al. 2008; Okamoto & Frenk 2009). The 
total number of massive dwarf galaxies with Umax > 30 kms“^ 
within the virial radius Uvir of the central galaxy is 6 in our Hydro 
simulation, which is almost half the value ( 11 ) of massive subha¬ 
los found in the DM0 simulation. Note that this corresponds to the 
mass range of the “massive failures” considered in Boylan-Kolchin 
et al. (2011, 2012). Still, our result is slightly higher than the to¬ 
tal number (4) of massive satellites in the Milky Way, including 


LMC and SMC, which have Umax above 30 kms“^ (Penarrubia 
et al. 2008a). Moreover, detailed variations from one main galaxy 
to another could, in principle, resolve the residual discrepancy. 

The sharp contrast in the number of dwarf galaxies between 
the DM0 and Hydro simulations highlights the critical role of bary- 
onic physics in galaxy formation, and it points to a potentially vi¬ 
able solution of the “missing satellite” and the “too big to fail” 
problems, in agreement with suggestions by some previous stud¬ 
ies (e.g.. Brooks et al. 2013; Sawala et al. 2014a; Mollitor et al. 
2015). 

3.2 Mass Functions of Subhalos 

In order to investigate effects of baryons on the subhalo mass, 
in Figure 6 we compare the subhalo mass M^ub (top panel) 
and the maximum circular velocity Umax (bottom panel) al z — 
0 of the matched pairs between the two simulations. As the 
fitting curve (black solid line) is below the diagonal dashed 
line (Msub(Hydro) = Maub(DMO), or Umax(Hydro) = 
Umax(DM0)), it is clearly seen that the majority of subhalos in 
the Hydro simulation are less massive than their counterparts in the 
DM0 simulation, similar to the subhalo abundance findings in Sec¬ 
tion 3.1. The subhalo mass function of the Hydro simulation peaks 
at ~ 5 X 10® /i“^M 0 , which is about a factor of 2 lower than the 
peak of the DM0 subhalo mass function at ~ 10^ /i“^M 0 . 

In the Hydro simulation, only massive subhalos can form 
stars. The minimum mass for subhalos to host star formation is 
log(Msub) = 7.5 (or Umax = 10 kms“^), although it may be 
affected by the resolution of the simulation. In the mass range 
between 10 ® /i“^M 0 and 10 ® /i“^M 0 where we have sufficient 
mass and spatial resolution, there is a mixture of “dark” sub¬ 
halos and bright satellites (subhalos that contain stars). Such a 
co-existence of dark subhalos and bright satellites implies that a 
linear Mhaio — AT* correlation, as commonly assumed in semi- 
analytical galaxy models and abundance matching techniques (e.g. 
Guo et al. 2010; Moster et al. 2013), may not hold in the dwarf 
galaxy regimes, since some massive halos do not host galaxies 
with stars. This would complicate the application of the abundance 
matching to dwarf galaxies (Garrison-Kimmel et al. 2014b; Guo & 
White 2014) and the assignment of galaxies to dark matter halos in 
A—body simulations. 

Another important parameter is the peak mass of each sub¬ 
halo, ATpeak, defined as the maximum mass attained by the pro¬ 
genitor before it was accreted by its host. Using the peak mass 
is currently the standard method in abundance matching or semi- 
analytical modeling when dealing with subhalos (e.g. Guo & White 
2014; Garrison-Kimmel et al. 2014b), since this quantity represents 
a physical state unmodified by the subsequent interaction between 
the subhalo and the host. Figure 7 shows a comparison of ATpeak 
from the DM0 and Hydro simulations. We find that subhalos be¬ 
low 10® h~^MQ generally have lower ATpeak in the Hydro simula¬ 
tion than in the DM0 simulation, and that subhalos with peak mass 
higher than 10® /i“^M 0 are able to form stars. However, there are a 
few “outliers”: two subhalos with peak mass above 10 ® /i“^M 0 re¬ 
main completely dark, while three subhalos with peak mass below 
lO®/i“^M 0 actually contain stars. The fitting of the data shows 
that, 

log ATpeak [Hydro] = l.lOlog ATpeak[DMO] - 1.05, (2) 

which means that the subhalo peak mass is ~30% lower in the 
Hydro simulation than in the DM0 simulation for ATpeak ~ 
10® /i“^M 0 , and ~44% lower for ATpeak ~ 10® /i“^M 0 . 
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log(V„^) [DMO] [km/s] 

Figure 6. Compaiison of subhalo properties of the matched pairs at 2 : = 0 
in both the DMO and Hydro simulations. The top panel shows the sub¬ 
halo mass Maul, distribution function, while the bottom panel shows the 
maximum circular velocity Dmax distribution function. The grey dots rep¬ 
resent all dark subhalos (subhalos that did not form stars), while the filled 
red circles represent bright satellites (subhalos that contain stars). The solid 
black line is the fitting curve for all the subhalos (including both dark 
and star-forming ones), while the dashed line indicates Mgub (Hydro) = 
Msub(DMO) in the top panel, and ituiax (Hydro) = i;uiax(DMO) in the 
bottom panel. 


Our results suggest that the impact of hydrodynamics on the 
halo mass could he easily amplified in the early growth stages when 
the halos increase their mass exponentially (Bosch et al. 2014; Cor¬ 
rea et al. 2015), and that the assumption of a monotonic relation be¬ 
tween stellar mass of a galaxy and its peak mass in the abundance 
matching technique is not valid. Hydrodynamic simulations should 
hence be employed for a more reliable study of the properties of 
dwarf galaxies, also as suggested by previous work (e.g. Sawala 
et al. 2014a; Velliscig et al. 2014). 



Figure 7. A comparison of the subhalo peak mass Mpeak (th^ maxi¬ 
mum mass attained by the progenitor before it was accreted by its host) 
of matched pairs in both the DMO and Hydro simulations. The black filled 
circles are the “dark” subhalos and the red filled circles represent bright 
satellites. The black solid line is the fitting curve of all subhalos, while the 
dashed line indicates Mp^ak (Hydro) = Mpeak (DMO). 

3.3 Dark Matter Distribution and Density Profile of the 
Main Host 

To study the effects of baryons on the DM distribution of the main 
host, we compare the DM shape and density profile of the main 
host in both simulations. We apply a principal component analy¬ 
sis to the DM halo and compute the three axis parameters a, b and 
c based on the eigenvalues of the moment of inertia tensor for all 
the DM particles within a given shell (following the method by 
Zemp et al. 2011). The halo shape can be quantified by the in¬ 
termediate to major axis ratio, b/a, and the minor to major axis 
ratio, c/a, as shown in Figure 8 (top panel), which compares the 
axis ratios at different distances from the galactic center in both the 
DMO and Hydro simulations. The shape of the DM halo differs sig¬ 
nificantly between the two simulations. In the inner region within 
10 /i“^kpc, it is triaxial with the triaxiality parameter (defined as 
T — [a^ — b^]/[a^ — (?\ as in Zemp et al. 2011) T ~ 0.9 in the 
DMO simulation, but in the Hydro simulation it is close to an oblate 
spheroid with 6/a ~ 1 and c/a > 0.7. The difference in halo shape 
between the DMO and Hydro simulations continues towards larger 
galactic distance out to ~ 100/i~^kpc, but they converge at the 
virial radius r ~ 200 6,“^kpc. These results are in good agreement 
with previous studies (Springel et al. 2004; D’Onghia et al. 2010a; 
Vera-Ciro et al. 2011; Zemp et al. 2011; Bryan et al. 2013), and 
demonstrate that the impact of baryons on the dark matter halo 
shape is significant up to the virial radius, a spatial scale much 
larger than that of the stellar disk or the central galaxy. We note that 
a similar effect is seen in gas-rich mergers of galaxies, where gas 
inflows (Barnes & Hemquist 1991) and nuclear starbursts (Mihos 
& Hemquist 1996) can significantly modify the shapes and orbital 
properties of remnants (Barnes & Hemquist 1996). 

In addition to the halo shape, remarkable differences between 
the Hydro and DMO simulations are also present in the DM den¬ 
sity profile of the main host. Adiabatic contraction of the DM dis¬ 
tribution by baryonic processes is a well known effect. Historically, 
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Figure 8. Comparison of the shape of the main dark matter halo from the 
DM0 (in black) and Hydro (in red) simulations. Top panel: the intermediate 
to major axis ratio b/a (solid lines) and the minor to major axis ratio c/a 
(dashed lines) as a function of distance to the galactic center. Bottom panel: 
the triaxiality T, defined as T = {a^ — b^)/{a‘^ — c^), of the main dark 
matter halo as a function of galactic radius. The dotted lines indicate T = 
0.25 and T = 0.75, marking the transitions from oblate/prolate to triaxial 
halo shapes. The plot shows that in the inner region the DM halo is more 
spheroidal in the Hydro simulation, but it is triaxial in the DM0 simulation. 

adiabatic contraction effects were calculated analytically based on 
the assumption of circular orbits and conservation of angular mo¬ 
mentum. Gnedin et al. (2004) refined the calculation by consider¬ 
ing the eccentricities of the orhits in a more realistic cosmological 
context. Even with such a modification, the density profile of DM 
in the inner region could be overestimated without the inclusion 
of more baryon physics other than cooling and star formation. Re¬ 
cently, Marinacci et al. (2014a) showed that an enhancement of the 
DM density in the inner region was present for most of the eight 
Milky Way-size halos studied in their simulations (however at a 
lower resolution than that used here). 

To illustrate the effect of adiabatic contraction on the DM 
distribution, we present in Figure 9 the spherically-averaged DM 
density profile of the main halo both from the DM0 and Hydro 
simulations (top panel), as well as a comparison with the adiabatic 
contraction calculation (bottom panel) based on the DM0 density 
profile using the contra^ code (Gnedin et al. 2004, 2011). This 
code calculates the response of a DM distribution to the condensa¬ 
tion of baryons. We input the DM distribution at 2 = 0 from the 
DM0 simulation as the state before contraction, and the baryonic 
mass distribution at 2 = 0 from the Hydro simulation as the source 
of adiabatic contraction. Since the contraction of DM is naturally 
followed in the Hydro simulation, we can compare the DM density 
profile from the Hydro simulation with the theoretical expectation 
from CONTRA. The lower panel of Figure 9 shows the expected en¬ 
hancement of DM from the DMO simulation (dash-dotted line) as 
if it would host the same galaxy produced by the Hydro simulation. 
The red solid symbols show the DM density profile measured in 
the Hydro simulation.This plot shows a significant enhancement of 

^ http;//dept.astro.Isa.umich.edu/-ognedin/contra. 



r/r 



r/r 

Figure 9. Comparison of the radial dark matter density profile of the main 
halo in our different simulations. Top panel: The halo DM density profiles 
from the DMO (black solid curve) and Hydro (red solid curve) simulations, 
respectively. A factor of is applied to the density profile from the 

DMO simulation. Bottom panel: The enhancement of DM from the Hydro 
simulation (red solid curve) in comparison with the adiabatic contraction 
calculation (blue dashed-dotted curve) based on the DMO density profile 
using the CONTRA code (Gnedin et al. 2011). The error bars are based on 
the Poisson en'or \/N, where N is the number of DM particles in each 
radial bin. The vertical dotted line in both plots marks the position at which 
the gravity force equals its exact Newtonian form, r = 2.8 e, where e is the 
gravitational softening length. An enhancement of the DM concentration 
due to adiabatic contraction is present in the Hydro simulation in the inner 
region up to r ~ O.lrviy. 


the DM in the inner region out to r ~ 0.1 Pyir in the Hydro simu¬ 
lation, which is consistent with the expected adiabatic contraction. 
Interestingly, in the very central region, we see more DM enhance¬ 
ment in the Hydro simulation compared to the result from contra. 
While it is possible that the DM distribution in the Hydro simu¬ 
lation follows a much more complex evolution than the simplified 
analytical calculation in contra, we also note that the “bump” in 
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the lower panel of Figure 9 may be subject to numerical uncer¬ 
tainty as it occurs on a scale very close to the spatial resolution of 
the simulations. 

In our Hydro simulation, the response of DM due to gas cool¬ 
ing and condensation is consistent with the expectation of adiabatic 
contraction. As demonstrated in Marinacci et al. (2014a), this holds 
true for the majority of the MW-size halos in the hydrodynamic 
simulations. The agreement between our simulation and contra 
fitting suggests that the latter provides a good description of the DM 
distribution in L* galaxies. Since contra has been calibrated using 
a suite of different hydrodynamic simulations (Gnedin et al. 2011), 
the agreement between our result and contra therefore reflects a 
consistent response of dark matter with respect to gas cooling and 
condensation, despite substantial differences in how feedback was 
modeled. 

Interestingly, no dark matter core is formed in our simulated 
halo. We note that recent simulations of halos more massive than 
10^^ Mq by Maccio et al. (2012) and Mollitor et al. (2015) show 
flattened DM distributions within the inner 5 kpc. Compared to 
our simulations, these studies employed a different feedback model 
featuring a more bursty stellar feedback. While feedback processes 
such as energy and momentum from SNe explosions may hence 
change the contraction of DM in the central region, such effects 
sensitively depend also on how the feedback injection is modeled 
in detail. 


3.4 Dark Matter Distributions and Density Profiles of 
Subhalos 

Similar to the main halo, we identify subhalos in both the DM0 
and Hydro simulations using the AHF group finder. We first locate 
the center of mass from the output of AHF for each matched object, 
then compute the spherically averaged density for DM particles of 
each subhalo based on the particle locations in the original snap¬ 
shot. Figure 10 shows the DM density profiles of 6 matched sub¬ 
halos at 2 = 0 from both simulations, as well as the corresponding 
circular velocity curves. These subhalos cover a wide range of total 
mass, from massive bright dwarf galaxies to “dark” subhalos. The 
circular velocity of each subhalo is calculated as y/GM{< r)/r, 
where M{< r) is the enclosed mass within r. We further compute 
the contributions from DM and baryonic components to the rota¬ 
tion curve in order to determine whether the differences in Umax 
from the two simulations are due to highly concentrated baryons 
or a genuine response of DM to baryonic processes. For the DM0 
simulation, we assume the distribution of baryons follows that of 
the DM but differs in mass by a factor of Ob/Dm- We note, how¬ 
ever, that the most massive subhalo. Sub 10, contains more bary¬ 
onic mass in the Hydro simulation than in the DMO simulation, 
while the least massive subhalos (such as sub 162 and sub 244) 
contain much less baryonic mass in the Hydro simulation than ex¬ 
pected based on the DMO simulation. 

From Figure 10, the computed density profiles of the subha¬ 
los from the Hydro simulation match their counterparts from the 
DMO simulation quite closely. However, the local distribution of 
DM is better probed by the rotation curves as they depend sensi¬ 
tively on the mass enclosed within a certain radius. As shown in 
the right panels of Figure 10, significant differences are evident 
between the rotation curves from these two simulations, in partic¬ 
ular with respect to the contribution of DM as represented by the 
solid curves in the right panels. For the first 3 subhalos (Sub 10, 
19 and 43), the contribution from the DM in the Hydro simulation 
is higher than that in the DMO simulation. They contain slightly 












Figure 10. A comparison of the dark matter density profiles (left column) 
and the circular velocity curves (right column) of 6 matched subhalos from 
the DMO (in black) and Hydro (in red) simulations . As in Figure 9, the 
DMO density profiles are multiplied by a factor of , and the dashed 

vertical line indicates r = 2.8 e. For each subhalo, its ID, Umax and bary¬ 
onic mass Mb are listed. In the right panels, the total circular velocity, and 
contributions from the DM and baryons are represented by the dotted, solid 
and dashed curves, respectively. The distribution of baryons in the DMO 
simulation is assumed to follow that of the DM multiplied by a factor of 
Db/Dm- 
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more dark matter in the Hydro simulation than in the DMO one, 
showing some mild contraction. In addition, the amount of con¬ 
traction in these three subhalos varies with their total mass, with 
Sub 10 showing the strongest contraction and sub 43 the weakest. 
On the other hand, the other three subhalos, sub 90, 162 and 244, 
show slightly reduced DM concentrations in the Hydro simulation 
compared with the DMO one. 

In the inner region, we have not found clear signs of DM cores 
in these subhalos. However, the absence of DM cores could poten¬ 
tially be due to a limited mass and spatial resolution of the sim¬ 
ulations. The gravitational softening length in our simulations is a 
factor of 2 to 3 larger than in the most recent high resolution simula¬ 
tions focusing on (isolated) dwarf galaxies (Govemato et al. 2012; 
Teyssier et al. 2013; Madau et al. 2014; Onorbe et al. 2015). 


4 THE IMPACT OF BARYONS ON THE EVOLUTION OF 
SUBHALOS 

4.1 Importance of Individual Physical Processes 

Baryonic processes play a critical role in the formation and evo¬ 
lution of galaxies. In this study we focus on three major mecha¬ 
nisms related to baryons that impact the mass distribution in galax¬ 
ies: adiabatic contraction, reionization, and tidal disruption. Adia¬ 
batic contraction due to gas cooling and condensation leads to an 
increase of the density in the galaxy center. Reionization not only 
ionizes and evaporates gas from galaxies, but can also prevent gas 
accretion from the intergalactic medium (IGM). Tidal truncation 
from gravitational interactions can result in the removal and redis¬ 
tribution of both dark and baryonic matter components. 

The relative impact of these process on galaxy properties and 
their evolution depends on the galaxy mass. Based on the prop¬ 
erties found in Section 3, the subhalos in our simulations can be 
categorized into three main groups according to their nmax (or al¬ 
ternatively mass) at z = 0. The first group consists of massive 
subhalos with Umax > 35 kms“^ (Msub > 4 x 10® M©), where 
adiabatic contraction tends to increase the amount of DM within 
the virial radius. These subhalos usually form stars and are there¬ 
fore “bright”. The second group are the least massive ones with 
■Umax < 20 kms“^ (Msub < 10® Mq), and are mostly “dark” 
with little or no star formation. These small subhalos are affected by 
reionization. The third group are subhalos with intermediate masses 
with Vmax ~ 20 kms“^ - 35 kms“^. These subhalos show signs 
of a competition between adiabatic contraction and tidal dismp- 
tion. While adiabatic contraction is able to increase the Umax of 
these subhalos, they often suffer from strong tidal effects in the 
Hydro simulation that remove both DM and baryonic mass, thus 
effectively reducing Vmax once they are close enough to the cen¬ 
tral galaxy. In what follows we will show some examples of these 
evolutionary paths. 
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Figure 12. Top panel: Comparison of the radially enclosed DM mass pro¬ 
files of Sub 10 from different simulations (top panel). The red solid, black 
dash-double-dotted and blue dashed-dotted curves represent the Hydro sim¬ 
ulation, the DMO simulation and CONTRA calculation, respectively, and the 
total enclosed baryonic mass (gas and stars) is also shown (blue tilled cir¬ 
cles). A factor of is applied to the density profile from the DMO 

simulation. This plot shows that Sub 10 retains more mass in the outer re¬ 
gion in the Hydro simulation than its DMO counterpart. Bottom panel: The 
radially enclosed DM mass profiles of Sub 10 in the DMO simulation at 
different evolution times. This subhalo has a pericentric passage at redshift 
2 = 0.10 in both simulations. In the DMO simulation, this subhalo experi¬ 
enced substantial mass loss in the outer region which explains the difference 
between the two mass profiles in the top panel beyond 10 kpc. 


4.1.1 The Role of Adiabatic Contraction 

As we have seen from Figure 10, three subhalos (Sub 10, 19 and 43) 
show signs of contracted DM in the Hydro simulation. Subhalo 10 
is the most massive among them, with Umax = 73kms“^ (in the 
Hydro simulation). Figure 11 shows the assembly history and dy¬ 
namical evolution of Sub 10. Not surprisingly, this object is able to 
fuel star formation continuously from high redshift to the present 
day, as shown by the “blue” stars in the composite images of the 


figure. Both the Hydro and DMO simulations produce similar tra¬ 
jectories for Sub 10, which is simply a single fly-by at 2 ~ 0.3. 
A slight reduction of the total mass and Umax is evident in both 
simulations due to this fly-by; but overall, the effect of baryons on 
this object is dominated by adiabatic contraction, making it more 
resilient to disruption. It builds up mass steadily, and it reaches a 
total mass of 2.1 x 10^® at 2 = 0 in the Hydro simulation, 

about 50% more massive than its counterpart in the DMO simula- 
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Figure 11. Time evolution of the most massive subhalo (Sub 10) in the simulations. Top panels: projected density map of the stellar component from redshift 
2 : ~ 5 to 2 : = 0. The composite images have used RGB colors mapped from K, B and U bands, respectively. The lower panels show (from left to right) the 
orbit of Sub 10 in the x-y plane and its distance to the center of the MW halo, as well as the evolution of its total mass and maximum circular velocity Dmax- 
As in previous plots, the DM0 and Hydro simulations are represented by black and red colors, respectively. In the last two panels, a vertical line represents 
the end of reionization in the Hydro simulation. 
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tion. It is interesting to note that the Umax of Sub 10 in the Hydro 
simulation is substantially larger than that in the DMO simulation 
throughout the time despite that its total mass is almost identical in 
the two runs. This suggests that the contraction of DM within Sub 
10 has been established at an early stage if adiabatic contraction is 
indeed responsible for the enhanced DM density profile shown in 
the top panel of Figure 10. Moreover, in the Hydro simulation there 
is a non-negligible contribution of the baryons to the total gravita¬ 
tional potential of Sub 10 in the inner regions which also helps 
explaining the different values of Umax (see Figure 12). 

To demonstrate the effect of adiabatic contraction in subhalos. 
Figure 12 shows a comparison of the radially enclosed DM mass of 
Sub 10 from both Hydro and DMO simulations against prediction 
from the contra code (top panel), and its evolution from the DMO 
simulation (bottom panel). In the inner region between 1 and 4 kpc, 
the Hydro simulation and the contra calculation show similar en¬ 
hancements of the enclosed DM mass by ~ 25% due to baryons. 
Beyond 4 kpc, the enclosed DM mass in the Hydro simulation is 
consistently higher than the expectation from contra. 

We find that the discrepancy between the Hydro and DMO 
simulations is mostly due to tidal removal of the loosely bound ma¬ 
terial in the outer parts of Sub 10 in the DMO simulation, as evi¬ 
denced by the change of radial enclosed DM mass with time in Fig¬ 
ure 12 (bottom panel). Sub 10 experiences continuous mass loss in 
the outer region after the pericentric passage at redshift 2 = 0.10. 
Although mass loss also occurs in Sub 10 in the Hydro simula¬ 
tion, as shown in the mass evolution in Figure 11, the amount is 
much smaller than that in the DMO simulation. These results show 
that when baryonic effects are included, massive systems similar 
to Sub 10 become more resilient to tidal disruption since adiabatic 
contraction and the presence of baryons in the inner regions tend to 
increase the binding energy. 

Similar effects of adiabatic contraction could also help explain 
the survival of the bright satellites of galaxy clusters reported for 
the Illustris Simulation (Vogelsberger et al. 2014a). It was found 
that those galaxies (with stellar mass ~ 10^° much more 

massive than Sub 10, the largest dwarf in our simulation) are more 
resilient to tidal disruption in the central cluster regions than satel¬ 
lites in pure A—body simulations due to the increased concentra¬ 
tion of DM and stellar components, in agreement with our findings 
in this study. 

The role of adiabatic contraction becomes progressively less 
important for lower mass dark matter halos/subhalos, as we have 
shown in Figure 10. In particular, for subhalos with Umax < 
35 km the total amount of baryons (mostly in the form of cool 
gas and stars) no longer plays a substantial gravitational role in 
these systems. We also note that our simulations include only a 
few massive subhalos, so that some scatter in their properties is 
inevitable. A more accurate estimate of the galaxy mass at which 
adiabatic contraction turns ineffective will require a much larger 
sample of (massive) subhalos than the one analyzed in this study. 

4.1.2 The Role of Reionization 

The epoch of reionization is an important landmark event in cos¬ 
mic history during which photons from young stars or accreting 
black holes ionize the neutral hydrogen. The latest results of Planck 
(Planck Collaboration et al. 2015) indicate that the Universe was 
50% reionized at 2 « 9, while Gunn-Peterson absorption features 
in quasar spectra suggest that reionization began as early as 2 ~ 14 
and ended at 2 ~ 6 (Fan et al. 2006). 

The majority of the low-mass subhalos in our simulations are 


unable to form any stars due to reionization, thus staying dark. In 
Figure 13, we show the evolutionary histories of five of such dark 
subhalos. Sub 149, 252, 347, 277 and 382. The trajectories show 
that they have been recently accreted onto the main halo. For each 
subhalo, the trajectory from the Hydro simulation is close to that 
of the DMO simulation, albeit with some small deviations. This 
is expected since the gravitational potential of the central host is 
modified in the hydrodynamic simulation due to the presence of 
the stellar disk, and the gas ram pressure, which is absent in colli¬ 
sionless A-body simulations, introduces some additional offsets in 
orbital phase space. 

The growth history of each subhalo shows remarkable differ¬ 
ences between the DMO and Hydro simulations after 2 = 6. The 
subhalo mass from the Hydro simulation is consistently lower than 
that from the DMO simulation. It is clear that the gas content of 
the subhalos declines rapidly after reionization, because the shal¬ 
low gravitational potential of these objects both fails to retain the 
heated gas and to accrete new gas from the IGM. 

A comparison of the growth histories between Figure 13 and 
Figure 11 shows that only the more massive subhalos are able to re¬ 
tain their gas after reionization, likely due to the fact that the dens¬ 
est gas regions in these objects can still reach the critical density 
needed for self-shielding from the ionizing UV background. These 
subhalos also have a sufficiently deep gravitational potential to ac¬ 
crete new gas and sustain star formation. These results are consis¬ 
tent with those obtained by Onorbe et al. (2015) for simulations of 
field dwarf galaxies, and demonstrate that reionization significantly 
suppresses the formation of low-mass galaxies. 

4.1.3 The Role of Tidal Disruption 

The third group we consider consists of bright satellites 
(intermediate-mass subhalos with 20 km s ^ < '^max < 

35 km s“^) similar to the dwarf spheroidal galaxies near the MW. 
We have already shown in Figure 3 the effect of increased tidal dis¬ 
ruption of subhalos in the Hydro simulation, which resulted in a 
lower subhalo abundance in the inner region of the main halo. Be¬ 
low we will examine individual subhalos and how they are shaped 
by tidal forces. 

In Figure 14, we map the distribution of the DM components 
of two satellites. Sub 436 and 244, at 2 = 0, both in the Hydro 
and DMO simulations. For the Hydro simulation we also examine 
the stellar distribution. The original members of each subhalo are 
identified at the redshift when their Umax reached the peak value. 
Overall, the trajectories of the matched subhalos are similar in the 
Hydro and DMO simulations, but with some subtle differences. For 
example, the current positions of the two subhalos in the DMO sim¬ 
ulation lag slightly behind those in the Hydro simulation. This may 
be caused by the deeper potential well of the main halo from a more 
contracted DM density profile and the presence of the stellar disk, 
which accelerate the subhalos close to the host to move at a slightly 
higher speed in the Hydro simulation than in the DMO one. How¬ 
ever, the most striking feature visible in the figure is the presence 
of tidal debris and streams of both DM and stars. This is clear ev¬ 
idence of strong gravitational interactions and tidal truncation of 
these two satellites. 

Figure 15 shows the evolutionary paths and growth histories 
of four subhalos (Sub 220, 162, 244, 436) that are massive enough 
to retain gas to fuel star formation in a continuous manner after 
reionization. The typical stellar mass of these subhalos falls in the 
range of 10® — 10^ at 2 = 0, similar to that of the classi¬ 

cal dwarf spheroidal galaxies in the Local Group (Boylan-Kolchin 


(© 2015 RAS, MNRAS 000, 1-22 



Baryonic impact on dark matter 13 




















Figure 13. The evolution of five low-mass subhalos (‘Umax < 20 kms~^) from the simulations (namely Sub 149, 252, 347, 277 and 382) in terms of their 
orbits in the x-y planes and their distance to the center of the MW halo at different redshifts (left two columns). Also shown are their growth histories in 
mass and circular velocity as a function of redshift (right two columns). As in previous plots, the DM0 and Hydro simulations are represented by black and 
red colors, respectively. In the two columns on the right, the total, gas and stellar masses are represented by solid, green dashed-dotted, and blue dotted lines, 
respectively, and the vertical dashed line indicates the end of reionization at 2 : = 6. 
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Figure 14. Projected maps of the dark matter distribution of two bright satellites (intermediate-mass subhalos with 20kms“^ < t/max < 35kms“^), 
namely Sub 436 (top panels) and 244 (bottom panels), at z = 0 from both the Hydro and DMO simulations. The left two columns show projections in the x-y 
plane, while the right two columns show the projections in the x-z plane. The blue dots represent the dark matter particles, while red dots show stars in the 
Hydro simulation. The presence of both dark matter and stellar streams in these plots are clear signs of tidal truncation. 


et al. 2012). These dwarfs have a total mass well above 10® h~^MQ 
before infall to the main host, and they can continue to accrete gas 
from the IGM after z = 6. 

These subhalos experienced strong gravitational interactions 
with their main host before the final plunges, as shown in their tra¬ 
jectories. These encounters tidally remove both baryons and DM, 
resulting in a steady reduction of the subhalo mass. The interactions 
also trigger episodes of starbursts, as shown in the growth curve of 
the stellar mass, and lead to distinct step-wise fluctuations in the 
velocity curve. The i/max curves experience stronger reductions in 
the Hydro simulation than in the DMO one. These are characteris¬ 
tic features of tidal forces during pericenter passages of the central 
galaxy. These findings are consistent with idealized simulations of 
tidal disruption of dwarf galaxies by Penarrubia et al. (2008b) and 
Arraki et al. (2014). 

In addition, the dwarf galaxies end up gas poor at 2 = 0 in 
the simulation, as shown in the evolution curve of the gas mass in 
Figure 15. All four dwarfs have experienced an abrupt loss in gas 
mass during their first infall. Gas is loosely bound to the subha¬ 
los compared to their stellar and DM components. However, tidal 
forces alone could not completely remove the gas. In principle, tidal 
torques may well funnel some gas into the center of these dwarfs 
and compress it to form a dense component which is difficult to 
strip (Mayer et al. 2006). However, ram-pressure, on the other hand, 
is able to efficiently remove the gas from these subhalos when they 
pass through the hot halo gas of the central galaxy (Mayer et al. 
2006; Wadepuhl & Springel 2011; Gatto et al. 2013; Arraki et al. 
2014). Ram-pressure stripping can transform these dwarf galaxies 
into gas-poor systems, and it may also induce other effects. For ex¬ 
ample, it was suggested by Arraki et al. (2014) that sudden gas loss 
in a dwarf galaxy due to ram-pressure could cause its dark matter 
to expand adiabatically and thus reduce Umax- However, this effect 
is small (< 10% in Umax) compared to the much larger reduction 
in Umax caused by tidal truncation. 


The trends of total mass evolution of subhalos in Figures 15 
show clear tidal disruption caused by the central galaxy. The sig¬ 
nature of tidal disruption is also clearly seen in the sharp and dis¬ 
tinct decrease of Umax- We show in Figure 16 the evolution of three 
subhalos with even shorter pericentric distances, comparable to the 
stellar disk size of ~ 25 kpc of the central galaxy (Marinacci et al. 
2014a). Not surprisingly, these subhalos experience even more sub¬ 
stantial reductions in both mass and Umax during their close en¬ 
counters with the host galaxy. In particular. Subhalos 3309 and 745 
show the largest reduction in Umax (~ 50% of their DMO val¬ 
ues) since they have passed the galactic center at much smaller dis¬ 
tances (~ 20 kpc) than the other subhalos. The mass loss in the 
stellar component is also higher for these three subhalos than those 
in Figure 15. 

The enhanced tidal disruption rate in the Hydro simulation 
is likely a combination of several gravitational effects, such as 
halo shocking from a rapidly varying potential which induces tidal 
shocks when the objects are on highly eccentric orbits, and disk 
shocking when they are passing in the vicinity of the stellar disk. 
In comparison, tidal stripping of material is a gentler process that 
does not increase the kinetic energy within the subhalo, at least 
when strong resonances are not operating (D’Onghia et al. 2009, 
2010b). To investigate the impact of tidal shocks on satellites, we 
follow the evolution of /r/umax, a ratio between the DM velocity 
dispersion a and the maximum circular velocity Umax a proxy of 
energy. In Figure 17 we show Subhalo 244 as an example during 
its infall journey into the main galaxy. Indeed, the internal energy 
of the subhalo, as indicated by a/umax, increases sharply when it 
passes the pericenter of its trajectory, demonstrating heating from 
tidal shocks during the close encounter. We have confirmed that the 
peaks of /r/umax are caused by the increase of a when the subhalo 
passes pericenter. We note, however, that only 64 snapshots were 
stored for the entire simulation, with the consequence that the time¬ 
sampling is not ideal for probing the subhalo trajectory at the time 
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Figure 15. The evolution of four bright satellites (intermediate-mass subhalos with 20kms“^ < Umax < 35kms“^) from the simulations (namely Sub 
220, 162, 244, and 436) in terms of their orbits in the x-y plane and their distance to the center of the MW halo at different redshifts (left two columns). 
Their growth histories in mass and circular velocity as a function of redshift are also shown (right two columns). As in previous plots, the DM0 and Hydro 
simulations are represented by black and red colors, respectively. In the two columns on the right, the total, gas and stellar masses are represented by solid, 
green dashed-dotted, and blue dotted lines, respectively, while the vertical dashed line indicates the end of reionization at z = 6. 


resolution required to clearly disentangle tidal from disk shocking 
or halo shocking. The discrete time sampling may also likely over¬ 
estimate the “minimum distance” plotted in these figures while the 
true “minimum distance” is attained between two snapshots. 

4.1.4 The Origin and Evolution of Bright and Dark Subhalos 

To understand the origin of “bright” (with stars) and “dark” (with¬ 
out stars) subhalos in the mass range of 10® —10® ^ Mq , we track 

their assembly histories in our simulations. In Figure 18, we show 


the distance of each subhalo to the central galaxy during its infall. 
Interestingly, “bright” and “dark” subhalos have different accretion 
paths. On average, the bright satellites are accreted into the host 
at an earlier redshift than the majority of the dark subhalos, and 
they typically undergo multiple passages through the main halo. 
Moreover, we follow the evolution of their Umax value and find that 
bright and dark subhalos have different trends as well, as shown in 
Figure 19. The dark subhalos experience a sharp decline in Umax (or 
mass) shortly after the end of reionization, while the bright satel¬ 
lites do not show such an immediate and dramatic suppression by 
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Figure 16. Same as Figure 15, but for three subhalos (namely Sub 608, 3309 and 745) with closest encounters with the central galaxy at distances comparable 
to the stellar disk size of ~ 25 kpc. As in previous plots, their orbits in the x-y plane and their distance to the center of the MW halo at different redshifts are 
shown in the left two columns, and their growth histories in mass and circular velocity as a function of redshift are shown in the right two columns. The DM0 
and Hydro simulations are indicated by black and red colors, the total, gas and stellar masses are represented by solid, green dashed-dotted, and blue dotted 
lines, respectively. Note that the redshift range log(l -|- z) is slightly narrowed compared to Figure 15 to highlight the temporal evolution of Umax, distance, 
and mass. 


reionization. In fact, most of them are able to retain the existing gas 
and replenish some of it long after 2 = 6, thus boosting their mass 
growth. In both cases, there is a reduction of Umax by ~ 17% on 
average at 2 = 0 in the Hydro simulation compared to the DMO 
one, highlighting the effects of baryonic processes on mass reduc¬ 
tion discussed in the previous sections. 

As demonstrated in Figure 19, reionization plays an important 
role in the formation of bright and dark satellites. The impact of 
reionization on these subhalos is mainly to suppress fresh gas ac¬ 
cretion from the IGM, as evidenced by a dip in the curve of the dark 
subhalos around log(l + 2 ) ~ 0.65 (2 ~ 3.5). However, a num¬ 
ber of the bright satellites in our simulation gain substantial mass 
even after 2 ~ 3 (see also Figure 13). Ricotti (2009) considered a 
scenario in which the low-mass halos (Vmax < 20kms“^) in the 
outer region of the MW are able to accrete low density IGM gas 
after 2 = 3 and form stars, once the mean temperature of the low- 
density IGM starts to decrease due to Hubble expansion, similar to 
what we find here. 

Our simulations show that bright satellites have a different ori¬ 
gin and evolutionary path from the dark ones. At early times, the 


bright satellites survive better than dark subhalos from reionization 
and they are able to retain gas and form stars afterwards. Moreover, 
they have an earlier infall time into the main galaxy than the dark 
ones, as also reported by other work (Sawala et al. 2014b). Our re¬ 
sults suggest that bright satellites may be biased tracers of the total 
subhalo population, as implied by observations of faint dwarfs of 
the Local Group (Weisz et al. 2015). 


5 DISCUSSION 

Thanks to significant progress in numerical modeling of galaxy for¬ 
mation and evolution, recent hydrodynamic simulations are becom¬ 
ing increasingly successful and are now able to reproduce many of 
the observed properties of galaxies in a more self-consistent man¬ 
ner. However, the complexity of the physical processes and the nu¬ 
merical methods inevitably entail uncertainties in the modeling. In 
what follows, we will compare our simulations with previous work 
and discuss the limitations of our model. 
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Figure 17. Effect of tidal shocks on the evolution of Subhalo 244. The black 
solid curves represent the evolution of cr/umax while the gray solid line 
represents the distance of Sub 244 to the center of the MW halo at differ¬ 
ent redshifts (in comoving units). Each peak of fr/wmax corresponds to a 
sharp increase of random kinetic energy due to tidal shock heating when 
the satellite passes through the pericenter. A blue vertical line indicates the 
infall time of this object at 2 : = 1.1, when it has first became a subhalo of 
the central host. 


5.1 Dark matter “cores” 

First of all, it is useful to compare the impact of baryons on the dis¬ 
tribution of DM in a MW-size halo found in our study with other 
works. The enhanced DM concentration in the inner region of the 
main halo seen in our simulation is consistent with adiabatic con¬ 
traction. Moreover, we do not find cored dark matter distributions 
on dwarf galaxy scales in our simulation. Recently, both Maccio 
et al. (2012) and Mollitor et al. (2015) have reported flattened DM 
distributions within 5 kpc from the galactic center in their hydrody¬ 
namic simulations of halos with similar mass. This difference can¬ 
not be attributed to insufficient numerical resolution as our mass 
resolution and gravitational softening length, which is finer than 
in Maccio et al. (2012), should give reliable results on kpc scales 
(Power et al. 2003; Springel et al. 2008). It is also unlikely to arise 
from differences in the hydrodynamics or gravity solvers, because 
grid-based codes such as RAMSES used in Mollitor et al. (2015) 
do not have the problem of intrinsic noise in SPH methods (Bauer 
& Springel 2012; Zhu et al. 2015) and use independent and differ¬ 
ent numerical methods than employed in Maccio et al. (2012). The 
most plausible cause for the difference lies in the feedback models. 
Here all three simulations have used a similarly large fraction of 
supernovae energy to drive outflows. However, our outflow model 
is less bursty than those of Maccio et al. (2012) and Mollitor et al. 
(2015), and this has been suggested as an important factor for mak¬ 
ing cores. We note however that Garrison-Kimmel et al. (2013) ar¬ 
gued that repeated blowout of gas is not necessarily more effective 
than a single blowout in reducing central dark matter densities. 

The outflow model used here is phenomenological and ties the 
wind launching velocity directly to the properties of subhalos. Once 
flagged as a wind, outflowing particles are temporally decoupled 
from hydrodynamics to prevent a disruption of the sub-resolution 
ISM phase. Thus, the small-scale creation of the wind in a Star- 


Figure 18. The distances (in comoving units) to the central galaxy of all 
subhalos in the mass range of 10® — 10® Mq at different times during their 
infall to the host, taken from the Hydro simulation. The red and black solid 
curves represent “bright” (with stars) and “dark” (without stars) subhalos, 
respectively, while the blue dashed curve indicates the virial radius of the 
main galaxy at different redshift. 


forming region and its interaction with the ISM is not followed in 
detail. Dalla Vecchia & Schaye (2008) have argued that such a de¬ 
coupled wind is less efficient than a coupled wind in driving strong 
turbulent motions on the scale of dwarf galaxies. A decoupled wind 
scheme could thus in principle miss some of the physical processes 
needed to generate DM cores, provided random bulk motions of 
gas are indeed able to produce a strongly fluctuating gravitational 
potential, as argued by Mashchenko et al. (2006); Pontzen & Gov- 
ernato (2012, 2014). However, Sawala et al. (2014a) reported that 
their simulated galaxies do not contain cores, even though they use 
a coupled wind model. This indicates that the decoupling feature of 
the wind model we used is not responsible for the absence of DM 
cores in our run. It thus remains interesting to see on what time 
scales the gravitational potential needs to fluctuate to generate DM 
cores. It is also possible that the delayed cooling mechanism used 
by Maccio et al. (2012) and Mollitor et al. (2015) overestimates the 
effect of bursty supernovae explosions (Agertz et al. 2013). 

Despite this difference between our simulation and others that 
predict DM cores, our results regarding the “dark” subhalos and the 
least luminous dwarfs should not be affected by the wind model, 
because these objects have the lowest star formation activities and 
hence an efficient removal of DM is energetically difficult (e.g. 
Govemato et al. 2012; Garrison-Kimmel et al. 2013; Madau et al. 
2014; Di Cintio et al. 2014; Ofiorbe et al. 2015). Indeed, Di Cin- 
tio et al. (2014) predicted that the most cored density distribu¬ 
tion is likely to be found in large halos with Umax = 50kms“^, 
and that DM profiles remain cuspy for dwarf galaxies with the 
least massive stellar population, in line with some recent hydro- 
dynamic simulations of field dwarf galaxies (Madau et al. 2014; 
Ofiorbe et al. 2015). The reduction on Umax at the low-mass end 
(Umax < 20kms“^, mostly “dark”) of the cumulative subhalo 
mass function in our simulation is in good agreement both with 
other SPH (Zolotov et al. 2012; Sawala et al. 2014a) and AMR 
(Mollitor et al. 2015) simulations, and we believe that reionization 
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Figure 19. Evolution of the maximum circular velocity of both dai'k (top 
panel) and bright (bottom panel) subhalos in the mass range of 10® — 
10^ Mq. These subhalos are identified in the Hydro simulation. In order 
to compare with the DM0 simulation and to identify the effects of baryons, 
the ratio 'Uinax,Hydro/'^max,DM of matched subhalos in both simulations 
is used. The vertical dashed line represents the end of reionization at red- 
shift 2 : = 6. The redshift bin size is 0.05 due to the limited number of 
output snapshots available for the simulations. 
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Figure 20. The radial distribution of the reduction of Umax 
(Avmax — (Vmax,Hydro Vmax,DM)/Vmax,DM) of matched Sub- 
halos between Hydro and DM0 simulations. The open black circles 
represent the dark subhalos that contain no stars, while each of the filled 
red circles represent individual bright satellites with its size proportional 
to the subhalo mass. The black dashed and the blue lines are simple 
linear fittings to the dark subhalos, and bright subhalos in the mass 
ranges 10^ M© < Mstar < 10® M© and 10® M© < Mgtar < lO"^ Mq, 
respectively, while the grey horizontal line provides a visual guide of no cor¬ 
relation between Avmax and galactic distance r. Overall, there is no clear 
radial dependence of Avmax of the subhalos in our simulations. Note the 
weak Avmax — t relation of subhalos with 10® Mq < Mstar < 10® Mq 
is due to two outliers at r = 50/i“^kpc. Once these two outliers are 
removed from the sample, the resulted radial dependence is as weak as the 
“dark” subhalos and those with 10® Mq < Mstar < 10^ Mq. 


is the primary culprit for the suppressed gas accretion rate from the 
IGM and the reduced number density of low-mass subhalos. 


5.2 Tidal disruption 

It was argued by Tollerud et al. (2014) that tidal forces play 
a marginal role in resolving the “too big to fail” problem as 
there is a lack of a strong radial dependence in the r — Umax 
relation for the M31 satellite galaxies. However, we find no 
strong radial dependence in the reduction of Umax in terms 
of (iimax, Hydro - «max.DMo)/«max,DMO in the Hydro simu¬ 
lation, as shown in Figure 20. Overall, the distrihution of 

(ttmax, Hydro nniax.DMO ) /ttmax.DMO IS rather flat, close tO 0.2, 
as a function of distance r. In this Figure, we further divide the 
subhalos into several subgroups according to their stellar mass and 
fit the results with simple linear functions. Only the subhalo sam¬ 
ple with 10® Mq < Mstar < 10® Mq shows a clear radial de¬ 
pendence. However, this strong radial trend is dominated by two 
subhalos. If we remove the two outliers around 50/i“^kpc, the 
radial trend disappears accordingly. For subhalos with 10®Mq < 
Mstar < IO^Mq, which is in the same mass range as the dSphs 
in Tollerud et al. (2014), there is no clear radial dependence for the 
reduction of Wmax. Hence, we conclude that a lack of a strong radial 
dependence cannot refute the role of tidal disruption. 

Tollerud et al. (2014) attributed the reduction of the number of 
subhalos to supernova feedback. However, we find that this plays 


a minor role in our simulation, especially for the subhalos with an 
intermediate Umax (see Fig. 14), in which the (maximum) circu¬ 
lar velocity is always larger in the Hydro case than in the DM0 
one before those objects are accreted into the central galaxy where 
stronger tidal disruption in the Hydro simulation comes into play. 
In our simulations, the tension between ACDM and observations 
of dwarf galaxies leading to the so-called “too big to fail” problem 
is largely alleviated by stronger tidal disruption, caused by an en¬ 
hanced DM concentration and the stellar disk of the central galaxy 
(See also Romano-Dlaz et al. 2010). The problem is further mit¬ 
igated by the fact that these subhalos are accreted by the host at 
a much earlier time than the average “dark” subhalos as shown in 
Figure 18. Thus, they have experienced a more extended tidal in¬ 
fluence from their host galaxy compared to the average subhalo. 
On the other hand, if the host halo in our simulation is slightly 
less massive, similar to the halo masses used by other groups (e.g. 
Guedes et al. 2011; Sawala et al. 2014a), one may not have a “too 
big to fail” problem at all to begin with (Wang et al. 2012). How¬ 
ever, it is unclear how tidal disruption or similar environmental ef¬ 
fects would work for the recently reported “too big to fail” problem 
of field dwarf galaxies (Garrison-Kimmel et al. 2014a; Papastergis 
et al. 2015; Klypin et al. 2014). 
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Figure 21. A comparison of our work with previous studies of the reduction 
of mass and and Umax. Top panel: Ratio of Mauij[Hydro]/Mguij[DMO] 
as a function of subhalo mass Mgui, [DMO] from our simulations and fitting 
relations from Sawala et al. (2013), Schaller et al. (2015) and Vogelsberger 
et al. (2014b). The filled symbols indicate our simulation data (red for bright 
satellites, black for dark subhalos), while the relations from Sawala et al. 
(2013) (GIMIC simulation) and Schaller et al. (2015) (EAGLE simulation) 
are shown as long-dashed and dashed-dotted lines, respectively. We include 
all the matched subhalos between the Illustris full physics run and lllustris- 
Dark without further separating them into subsamples as in Vogelsberger 
et al. (2014b, yellow line). Scatter from the fitting relations, shown as 1- 
(T error bars, is computed within each mass bin. Bottom panel: Ratio of 
Umax[Hydro]/Umax[DMO] as a function of UmaxlDMO]. The scatter of 
this relation, as indicated by l-tr error bars, is smaller than the top panel for 
Illustris. In both panels, a moving average of 200 data points is shown (black 
thick solid curve) to highlight the overall trend in the low-mass range. 


5.3 Mass reduction 

A number of previous works have examined the effects of baryons 
on halos in different mass ranges (e.g., Sawala et al. 2013; Schaller 
et al. 2015). In order to compare with these studies, we show in 
Figure 21 the reduction of mass and Umax from different simula¬ 


tions. The top panel compares the total mass reduction in the ha¬ 
los/subhalos as a function of subhalo mass from our simulations 
with the fitting relations from Sawala et al. (2013) (GIMIC sim¬ 
ulation) and Schaller et al. (2015) (EAGLE simulation), as well 
as all subhalos from the matched catalog between the Illustris 
and Illustris-Dark simulations (Nelson et al. 2015)^. To allow for 
a comparison with our results we extrapolated the fitting rela¬ 
tions of Sawala et al. (2013) and Schaller et al. (2015) down to 
10® Mq. Given the form of these relations this is equivalent to 
using a constant value of 0.65 and 0.73, respectively, for subha¬ 
los less massive than ~ 10® M©. Under this assumption, the plot 
shows that all studies are in good agreement for the mass reduc¬ 
tion of the central MW-size galaxy, as well as the low-mass end 
(< 5 X IO^Mq) where most of the dark subhalos and bright satel¬ 
lites in our simulations are located. There is substantial scatter in 
M 2 oo,Hydro/M 2 oo,DMO> which is both evident in the data points 
as well as the error bars of the Illustris results. We show a moving 
average of 200 data points with a thick solid black curve to high¬ 
light the overall trends in the low-mass range covered by the high 
resolution simulation used in this study. 

However, significant differences are present in the mass range 
between 5 x IO^Mq and lO^^M©. The mass ratio given by Sawala 
et al. (2013) or Schaller et al. (2015) is evidently lower than ours. 
Although we suffer from small-number statistics in our simulation, 
the mass excess in the Hydro simulation appears to be a real signal 
as also shown by a much larger galaxy sample drawn from the Il¬ 
lustris matched subhaloes. Considering that the universal baryonic 
mass fraction is ~ 16%, galaxies in Sawala et al. (2013) or Schaller 
et al. (2015) actually have less mass in the dark matter component 
in their hydrodynamic simulations. This led them to conclude that 
supernova feedback, which is strongly operating in this mass range, 
has also removed some DM, thus producing lighter halos. How¬ 
ever, we do not see DM mass reduction in this mass range in our 
simulation. As we have found in our earlier discussions on Sub 10 
(fmax = 73km s“^), we actually observe an increased dark matter 
mass in this regime due to adiabatic contraction. 

We note that the mass values returned by halo finders do have 
some uncertainties, which is evident in the large scatter of the raw 
data in the upper panel of Figure 21, while Umax is not as severely 
affected (e.g. Behroozi et al. 2015). In the lower panel of Fig¬ 
ure 21, we compare the ratio of Umax [Hydro]/Umax [DMO] from 
our simulations with that of the Illustris simulations. Both works 
show good agreement over the common Vmax. range, with subhalos 
above Umax ~ 35km (Maub > 5 X lO^M©) having higher 
Umax in the hydrodynamic simulation than in their DMO counter¬ 
part. Similar to the top panel, we show a moving average of 200 
data points at the low Umax end with a thick solid black curve to 
highlight the overall behavior, which is consistent with the Illustris 
relation. 

Since our simulations employ the same code and essentially 
the same physical models as the Illustris simulations, the differ¬ 
ences between our study and previous ones by Sawala et al. (2013) 
and Schaller et al. (2015) may owe to different implementation 
of feedback processes or different hydrodynamical methods used 
in these works. These comparisons thus highlight the need for 
comprehensive and systematic investigations similar to the Aquila 
Comparison Project, which is beyond the scope of this paper, but 
we plan to pursue it in future work. 


® Raw data from each simulation are available from http;//www. 
illustris-project.org/data/. 
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5.4 Implications for dark matter detections 

A^—body simulations are routinely used to model direct/indirect 
signals from Galactic dark matter structures. However, all the re¬ 
cent studies agree that baryons have significant impacts on dark 
matter distributions. The interplay between baryonic and DM dis¬ 
tributions has a direct impact on current attempts to indirectly mea¬ 
sure DM in the MW as well as in other galaxies. Traditionally, DM 
substructures are thought to be responsible for the observed radio 
flux-ratio anomaly in gravitational lensing (e.g. Mao & Schnei¬ 
der 1998). Recent studies by Xu et al. (2010, 2015) based on the 
A^—body simulations by Springel et al. (2008) and Gaoet al. (2012) 
concluded that DM substructures alone can not be the whole reason 
for radio flux-ratio anomalies. Xu et al. (2015) suggested that a sub¬ 
stantial improvement over the existing modeling of strong lensing 
would be to consider substructures in baryonic simulations that take 
the impact of baryonic processes on DM into account. Similarly, 
accurate predictions of the DM annihilation rate (The Fermi-LAT 
Collaboration et al. 2013) sensitively depend on the DM distribu¬ 
tion in the dwarf galaxies. The uncertainties in the DM distributions 
enter in factors describing the clumping of the DM density (propor¬ 
tional to the density squared) along the line of sight (see Geringer- 
Sameth et al. 2015). For example, Gomez-Vargas et al. (2013) have 
shown that the cross section estimated from 7 —ray photons from 
the galactic center can be affected by adiabatic contraction of DM 
by more than an order of magnitude. So far, these studies are based 
on inferences from high resolution A—body simulations. Hydrody- 
namical simulations change this picture and offer a more accurate 
and self-consistent physical modeling for analyzing direct and in¬ 
direct detection experiments. 


6 CONCLUSIONS 

The availability of hydrodynamical simulations allows us to re¬ 
assess some of the well-known potential issues of ACDM that 
have been identified with pure A-body simulations. In particu¬ 
lar, hydrodynamical simulations are able to self-consistently model 
the impact of baryonic processes on the DM distribution, a point 
which has previously often been ignored in galaxy formation stud¬ 
ies based on DM-only simulations. 

In this work, we analyze a high-resolution cosmological hy¬ 
drodynamic simulation on a moving mesh and compare it to its 
DM-only counterpart in order to study the DM distribution in a 
MW-like galaxy and its subhalos. This simulation uses essentially 
the same physical model employed in the Illustris simulation, and 
hence is consistent with a globally successful model for explain¬ 
ing the observed galaxy population. Moreover, the properties of the 
simulated galaxy are well converged with resolution, making our 
results numerically robust. Here we summarize our main findings: 

• We identify three physical processes induced by baryons 
that shape the overall distribution of DM in the main halo and 
its subhalos depending on mass: ( 1 ) adiabatic contraction due to 
gas cooling and condensation increases the DM concentration in 
the inner region of the main halo and in massive subhalos with 
fmax > 35 kms“^ (Msub > 4 X 10® h“^M©), making them 
more spherical in the inner regions; ( 2 ) reionization plays a crit¬ 
ical role in the formation and evolution of low-mass halos with 
Wmax < 20 kms“^ (Msuh < 10 ® h“^M©) by removing the gas 
from the halo and suppressing new gas accretion from the IGM, 
making them “dark” with little or no star formation; (3) strong 
tidal forces in the Hydro simulation effectively remove the stellar 


and DM components of intermediate-mass subhalos in the range of 
■Umax ~ 20 kms“^ - 35 kms“^ during their infall to the main 
galaxy, leaving behind tidal debris and streams of DM and stars. 

• As a result of these major effects from baryons, the total num¬ 
ber of subhalos in a MW-like galaxy and their total mass are sig¬ 
nificantly reduced in the hydrodynamic simulation compared to the 
DM-only one. Our results are in good agreement with observations 
of dwarf satellites in the MW galaxy, suggesting a viable solution 
to long-standing problems such as the “missing satellites” and “too 
big to fail” issues that arose from pure A-body simulations. 

• The ranking of subhalos based on either their peak mass or 
their present-day mass is modified in the Hydro simulation com¬ 
pared with the DMO run. A large fraction of subhalos with a peak 
mass below lO®h“^M 0 are found to host no stars. These findings 
suggest that the assumption of a monotonic relation between stellar 
mass and peak mass as commonly used in abundance matching is 
not strictly valid, and that effects of baryonic processes should be 
included in this modeling as well. 

Interestingly, our high-resolution hydrodynamic simulation 
does not produce cored DM distributions as observationally sug¬ 
gested in some low surface brightness galaxies, and unlike some 
other recent simulation models with very bursty feedback modes. 
More sophisticated hydrodynamic simulations are needed to fur¬ 
ther study this puzzling problem, and to improve the realism of 
the inclusion of feedback processes to determine whether they can 
indeed provide a solution to small-scale tensions identified in the 
ACDM cosmogony. 
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